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GRMHD Simulations of Magnetized Advection Dominated 
Accretion on a Non-Spinning Black Hole: Role of Outflows 
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ABSTRACT 

We present results from two long-duration GRMHD simulations of advection-dominated 
accretion around a non- spinning black hole. The first simulation was designed to avoid signif- 
icant accumulation of magnetic flux around the black hole. This simulation was run for a time 
of 200, 000GM/ c 3 and achieved inflow equilibrium out to a radius ~ 90GM/ c 2 . Even at this 
relatively large radius, the mass outflow rate M out is found to be only 60% of the net mass 
inflow rate Mbh into the black hole. The second simulation was designed to achieve substan- 
tial magnetic flux accumulation around the black hole in a magnetically arrested disc. This 
simulation was run for a shorter time of 100, OOOGM/c 3 . Nevertheless, because the mean 
radial velocity was several times larger than in the first simulation, it reached inflow equilib- 
rium out to a radius ~ 170GM/c 2 . Here, M out becomes equal to M B h at r - 160GM/c 2 . 
Since the mass outflow rates in the two simulations do not show robust convergence with time, 
it is likely that the true outflow rates are lower than our estimates. The effect of black hole 
spin on mass outflow remains to be explored. Neither simulation shows strong evidence for 
convection, though a complete analysis including the effect of magnetic fields is left for the 
future. 

Key words: galaxies: jets, accretion, accretion discs, black hole physics, convection, binaries: 
close, methods: numerical 



1 INTRODUCTION 

Black hole (BH) accretion occurs via at least two distinct modes : 
(1) a standard thin acc r etion disc dShakura & Sunvaevl 11973b 
iNovikov & Thornd Il973l : iFrank. King & Raind l2002n. and (2) 
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for reviews). Thin discs are 


present around stellar-mass 



and supermassive BHs that accrete at a substantial fraction ~ 
few-100% of the Eddington rate, while ADAFs are typically found 
at lower accretion rates M0 

The accreting gas in an ADAF is radiatively inefficient; hence 
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1 Actually, two distinct ADAF modes are possible, one in which optically 
thin two-temperature gas accretes with a highly sub-Eddington M, and a 
second in which very optically thick radiation-trapped gas accretes at rates 
well above the Eddington rate. We are concerned in this paper with the 
former kind of ADAF, which our GRMHD code is cap able of simulating. 
The latter variety of ADAF is referred to as a "slim disc" (Abramowic z et al.l 



an ADAF is also referred to as a radiatively inefficient accretion 
flow (RIAF). The low radiative efficiency, on top of the already 
low accretion rate, makes ADAFs highly underluminous and dif- 
ficult to observe. On the other hand, the vast majority of both 
stellar-mass and supermassive BHs in the universe are in the ADAF 
state, a notable example being Sgr A*, the supermassive BH at 
the center of our own Galaxy dNaravam Yi & Mahadevanl 1 19951 ; 
lYuan. Ouataert & Naravanll2003h . 

A simple one-di mensional model of gas dynamics in an ADAF 
(INaravan & Yil ll994) reveals two interesting complications. First, 
the Bernoulli parameter of the gas tends to be positive. This means 
that the gas is not gravitationally bound to the BH, or at best is 
only weakly bound. Therefore, an ADAF is likely to have power- 
ful jets and mass outflows, as recognized in the very first papers 
(Naravan & Yi 1994, 1995a). The connection between ADAFs and 
relativistic jets has become in creasingly clear over the years (e.g., 
INaravan & McClintockl 120081 and references therein). However, it 
is presently unknown whether or not ADAFs have quasi- or non- 



1988) and requires a radiation M HD code to simulate (see lOhsuga et al.l 
2009; Ohsuga & Mineshige 2011; and references therein). 
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relativistic winds, and if so how much mass they lose via these 

winds. 

Some authors (e.g.. iBlandford & Begelman1ll999l : iBegelmanl 

l2012h have suggested that winds in ADAFs are so powerful that 
the mass accretion rate Mbh on the BH is as much as ~ 5 or- 
ders of magnitude less than the mass supply rate M supp i y at the 
outer edge of the accretion flow, say at the Bondi radius. In effect, 
these authors took the Bernoulli argu ment for strong outflows pro - 
posed in the original ADAF papers dNaravan & Yilll994l Il995al) . 
and prostulated that ADAFs would have not just st rong outflows , 
but ov erwhelmingly strong outflows. Other a uthors (IOgilvielll999l : 
lAbramowicz. Lasota & Igumenshchevll2000h . however, argued that 
the Bernoulli parameter is not a good diagnostic for mass loss, es- 
peci ally in the case of viscous n o n- stead y flows. 

lYuan. Ouataert & Naravanl d2003h attempted to constrain 
the mass loss in Sgr A* using radio data on Faraday rota- 
tion (lAitken et al.ll2000l: lOuataert & Gruzinovll2000al : lAgolll2000l : 
iBower et"aD 120031 ; iMarrone et al.l 120071) . They concluded that, 
for this source, the decrease of M between the Bondi ra- 
dius and the BH is on the order of one to two orders 
of m a gnitude. More recently, a few studies (e.g ., lAllen et aD 
iMcNamarZ Rohanizadega n & Nulsenl 120111) have shown 



that many radio-loud active galactic nuclei require a power 
source comparable to or even greater than what Bondi accre- 
tion can supply. Even if the power source of the jet is BH 
spin energy, one still requires a signifi cant mass accretion rate 
on to the BH to tap this spin powe r dNaravan & Fabianll201ll : 
iTchekhovskov. Naravan & McKinnevI 1201 lh . Therefore, in the 
above radio sources, there cannot be significant mass loss between 
the Bondi radius and the BH horizon. 

The second potential complication in the dynamics of ADAFs 
is that the entropy gradient is la rge and highly unstab le according 
to the Schwarzschild criterion dNaravan & Yilll994l) . One might 
thus suspect that ADAFs will be very convective. On the other 
hand, the angular momentum profile has a stable gradient. It 
is thus not clear whether the flow is ultimately stable or un- 
stable to convection. An alytical models of convection-dominated 
accre t ion flows (CDAFs; LNarayan. Igumenshchev & Abramowiczl 
2000: IOuataert & Gruzino v 2000b) have bee n developed, but their 
relevance to real ADA Fs is unclear (see iNaravan et al.l 120021: 
iBalbus & Hawlevll2002l : for conflicting views). 

Both mass-loss and convection involve multi-dimensional 
flows, which are best studied via numerical simulations. In addi- 
tion, since the "viscosity" that dri ves accretion originates in th e 
magnetorotational instability (MRI. lBalbus & Hawlevli99lLll998h . 
magnetic fields play a critical role. This makes analytical studies 
even less tractable. Fortunately, multidimensional numerical MHD 
simulations are now feasible. Indeed, the limit of a non-radiative 
ADAF is relatively easy to simulate, since there is no radiation 
physics involved. Moreover, ADAFs are geometrically thick and 
are less demanding in terms of spatial resolution. We briefly review 
here the large literature on ADAF simulations. 

Early numerical simulations of ADAFs employed pseudo- 
Newtonian co des with purely hydrodynamic vis cosity. Pioneer- 
ing work by IStone. Pringle & Begelmanl dl999h indicated that 
such flows are convective and that a significant fraction of 
the inflowing mass near the equatorial plane flows out along 
the poles in a strong outflow. Similar results, viz., convection, 
equatorial inflow and bipolar outflow, w ere obtained also by 
llgumenshchev & Abramowiczl dl999Ll2000l) . In the latter paper, the 
authors found that bipolar outflows required high values of the vis- 
cosity parameter a, while low-viscosity models exhibited weaker 



outflow s but had strong convection. Very recently, I Yuan. Wu & Bui 
<2012bL see also lYuan. Bu & Wull2012ah have carried out 2D hy- 
drodynamic simulations of ADAFs which cover a very large range 
of radius and show fairly strong outflows. Most of the outflow- 
ing gas is bound to the BH in the sense that it has a negative 
Bernoulli parameter, yet it rea ches the outer boundary of the simu- 
lation without turning around. Ili. Ostriker & Sunvaevl(l2012l) have 
carried out hydrodynamic simulations of ADAFs including the ef- 
fects of bremsstrahlung cooling and electron thermal conduction. 

Although interesting, hydrodynamic a-viscosity simulations 
are ultimately not realistic since accretion flows have magnetic 
fields and MRI-driven turbulence. It is thus necessary to in- 
clude magnetic fields consistently. Pseudo-Newtonian magneto- 
hydrodynamic (MH D) simulations have been performed by a 
number of authors. Machida. Hayashi & Matsumotol d200Qh and 
iMachida. Matsumoto & Mineshige ~ l200ll) observed temporary 
outflows of mass in their MHD simulations and showed that sub- 
stantial accretion energy can be released in the vicinity of the BH 
via magnetic reconnection. They also claimed that the initial con- 
figuration of the magnetic field may play an important role in de- 
ter mining the mass outflow rate. Using axisymmetric (2D) mod- 
els, IStone & Pringlel d200ll) showed that significant outflows orig- 
inate at radii beyond r ~ 10 (we express length s in BH mass 
units: GM/c 2 ). Similarly. ISawlev & Balbusl (120021) observed out- 
flows for all radii outside the innermost stable circular orbit (ISCO), 
though they used a definition of inflow and outflow based on cy- 
clindrical coordinates (all other authors use spherical coordinates) 
which makes their outflow estimates somewhat ambiguous. 

Conv ective motions were e vident in MHD simulations per- 
formed bv lMachida et all d200ll) . indicating, according to the au- 
thors, that convection is a rather general phenome non in radiatively 
ineffic ient accretion flows. On the other hand, IStone & Pringlel 
d200ll) concluded that the turbulence seen in their MHD sim- 
ulations was driven by the MRI, not convection. Similarly, 
lHawlev & Balbusl d2002r) noted that, although their models were 
unstable according to the classical Hoiland criteria, the flows ap- 
peared not to be convective. On the other han d, a simulation by 
llgumenshchev. Naravan & Abramowiczl d2003l) . which was initial- 
ized with purely toroidal magnetic field, showed significant con- 
vection, and appeared to be similar to a CDAF. The same au- 
thors found that, if they initialized the simulation with a poloidal 
magnetic field, the disc structure was completely different from 
the toroidal case. The poloidal case led to a configuration in 
which the magnetic field strongly resisted the accreting gas, lead- 
ing to what the authors later called a "magnetically ar rested disc" 
(MAD, INaravan. Igumenshchev & Abramowiczl 120031). In a s eries 
of numerical MHD simulation s, Pen, Matzner "& Wond J2003h and 
IPang. Pen. Matzner. Green & L iebendorfer (2 01 lh found little evi- 
dence for either outflows or convection. Even though the entropy 
gradient was unstable the gas was apparently prevented from be- 
coming convective by the magnetic field. They coined the term 
"frustrated convection" to describe t his behavior. 

Beginning with the work of |De Villiers, Ha wlev & Krolikl 
(2003), accretion flows have been studied using general relativistic 
magneto-hydrodynamic (GRMHD) codes. iDe Villiers et all J2003I) 
observed two kinds of outflows: bipolar unbound jets and bound 
coronal flow. The coronal flow supplied gas and magnetic field to 
the coronal envelope, but apparently did not have sufficient energy 
to escape to infinity. The jets on the other hand were relativistic 
and escaped easily, though carrying very litt le mass. Jets have been 



studied in detail by a numbe r of author s (McKinnev & Gammie 
120041 : IDe Villiers etaD 120051 : iMcKinnevl l2006h . iBeckwith et al 
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(2008, 2003) and iMcKinnev & Blandford ( 2009) noted that the 
power emerging in the jets depended strongly on the assumed mag- 
netic field configuration. While dipolar fields produced strong jets, 
a qu adrupolar field led to only w eak, turbulent outflows. 

Tchekhovskov et al ]<En|) simulated a MAD system around 
a rapidly spinning BH, and obtained very powerful jets with en- 
ergy efficiency 77 > 100%, i.e., jet power greater than 100% of 
Mbhc 2 , where Mbh is the mass accretion rate on to the BH. Their 
work showed beyond doubt that at least some part of the jet power 
had to be extracted from the spin energy of the BH. The jet- spin 
co nnection for MAD systems has been explored in greater detail 
bv IMcKinnev. Tchekhovskov & Blandford ( l2012h . These authors 
coined the term "magnetically choked accretion flow" (MCAF) to 
describe the MAD configuration. 

Returning to the present paper, the goal here is to use GRMHD 
simulations of ADAFs around BHs to investigate the importance of 
mass outflows, and if possible convection. Our simulations are run 
for a longer duration than most previous work. The questions we 
address require us to analyze the properties of the accretion flow 
over as wide a range of radius as possible. The only way to obtain 
converged results over such large volumes is by running simula- 
tions for a very long time. We introduce a new measure of conver- 
gence, or more accurately a test of internal consistency. As per this 
criterion, our simulations give converged time- steady flows over a 
range of up to 100 in radius. This turns out to be still not as large as 
we would like. Nevertheless, it permits us to reach some interesting 
conclusions. 

Within the realm of ADAFs, we expect answers to depend on 
several factors. One important factor has already been mentioned, 
viz., the magnetic field topology in the accreting gas. The role of 
field topology for mass outflows (as distinct fro m relativistic jets) 
has be en largely unexplored. The recent work of IMcKinnev et all 
(120121) is one of the first studies in this area. 

In this paper we consider two distinct magnetic topologies 
and describe one long-duration simulation for each topology. In 
one simulation, we carefully arrange the initial seed magnetic field 
(which is later amplified via the MRI) such that the accreting gas 
does not become magnetically arrested despite the long duration 
of the simulation. We call this the ADAF/SANE simulation (where 
SANE stands for "standard and normal evolution"). In the second 
simulation, we set up the magnetic field topology such that the ac- 
cretion flow very quickly becomes magnetically arrested and then 
remains in this state for the duration of the run. We call this the 
ADAF/MAD simulation (where, as stated earlier, MAD stands for 
"magnetically arrested disc"). 

A second obvious parameter that will affect the properties 
of an ADAF is the spin of the central BH. Numerical studies of 
jets, for instan ce, clearly show that jet power correlates strongly 
with BH spin (IMcKinnev! 1 200d: I Tchekhovskov et alJl201 lL l2012l : 
iTchekhovskov & McKinnevll2012|). Observationally too, ther e is 
evidence for such a correlation iNaravan & McClintocklEoT2l) . In 
this paper we focus on the case of a non-spinning BH: a* = 
a/M — 0. We view such a system as the purest form of an ADAF, 
where the only available energy source is gravitational potential en- 
ergy released via accretion. By simulating an ADAF around a non- 
spinning BH using a GRMHD code, we can more easily relate our 
results to analytical studies as well as previous non-relativistic sim- 
ulations. In the future we plan to run long-duration GRMHD sim- 
ulations of ADAFs around spinning BHs. Those simulations will 
have two sources of energy, accretion and BH spin. By compar- 
ing them with the simulations described here we should be able to 
evaluate the role of BH spin. 



The plan of the paper is as follows. In §2, we briefly describe 
the simulation methods we employ, which are similar to those we 
have used in previous work. In §3, we discuss in detail our results 
from the ADAF/SANE and ADAF/MAD simulations, focusing in 
particular on mass outflows. In §4, we bring together the results of 
the previous sections and try to assess the nature of the accretion 
flow in the two simulations. In §5, we conclude with a discussion. 



2 DETAILS OF THE SIMULATIONS 

2.1 Computation Method 

The simulati ons described here were done with the 3D GRMHD 
code HARM dGammie McKinnev & Tothll 20031 : lMcKinne\ll20()3 : 
IMcKinnev & Blandfordl 20091) . which solves the ideal MHD equa- 
tions of motion of magnetized gas in the fixed general relativistic 
metric of a stationary BH. The equation of state of the gas is taken 
to be u = p/(T — 1), where u and p are the internal energy and 
pressure, and T is the adiabatic index. The code conserves energy 
to machine precision, hence any energy lost at the grid scale, e.g., 
through turbulent dissipation or numerical reconnection, is returned 
as entropy of the gas. There is no radiative cooling. The code works 
in dimensionless units where GM = c = 1. Thus, all lengths and 
times in this paper are given in units of GM/c 2 and GM/c 3 , re- 
spectively. 

A key feature of our simulations is the extremely long run 
time: 200, 000 time units for the ADAF/SANE simulation, and 
100, 000 time units for the ADAF/MAD simulation. To avoid spu- 
rious signals reaching the region of interest from the boundary of 
the simulation, our grid extends out to a very large radius ~ 10 5 . 
At the same time, we require good resolution in the inner regions 
in order to study the structure of the flow. To satisfy both require- 
ments, we use a grid with 256 cells in the radial direction, where the 
cells are distributed uniformly in log r at smaller radii and spaced 
hyper-logarithmically near the outermost radii. 

In the direction, we employ 128 cells, distributed non- 
uniformly so as to provide adequate resolution both in the geomet- 
rically thick equatorial region, where the bulk of the gas accretes, 
and the polar region, where a relativistic jet might flow out. In order 
to follow suc h a jet as it collimates at lam e distance, we use the grid 
developed bv T chekhovskov et al J (120111) in which the resolution 
near the pole increases with increasing radius (see Fig. 

EE 

Finally, we use a uniform grid of 64 cells in the azimuthal 
direction, covering the full range of <j) from to 2tt. 

2.2 Initial Conditions 

The fluid initially rotates around the BH in a torus 

in hydrostatic equilibrium: a ^Polish doughnut" 

(iKozlowski Jaroszvnski & Abramowiczll 19781) . The ADAF/SANE 
and ADAF/MAD simulations begin with the same torus. It has 
inner edge at n n = 10 and extends to r ~ 1000 (Figs. [2 [3]). The 
angular momentum of the torus is constant inside r break = 42. 
Outside rbreak, the angular momentum is 71% of the Keplerian 
value and is constant on von Zeipel cylinders. The entropy is 
constant everywhere, p/p r = 0.00766, and the Bernoulli is small 



2 As it happens there is no significant jet in the simulations described here. 
However, we plan to use the same grid setup and initial conditions in future 
work with spinning BHs, where we do expect to see strong jets. 
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Figure 1. Poloidal plane of the grid used in the simulations, shown at two zoom levels. 




Figure 2. Initial configuration of the ADAF/SANE simulation. The top two panels show the mid-plane density and the magnetic flux threading the equatorial 
plane as a function of radius. Note the extended size of the initial torus, which is required for the extremely long duration of this simulation. Note also the 
multiple oscillations in the magnetic flux, which prevents the accreting gas from reaching the MAD state. The lower two panels show the logarithm of the 
density p and the gas-to-magnetic pressure ratio (3 of the initial torus in the poloidal plane. 
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Figure 3. Similar to Fig.[2]but for the ADAF/MAD simulation. The main difference is that the torus here has a single loop of field centered at radius r = 300. 
As a result, accretion causes magnetic flux of one sign to accumulate around the BH, leading to the MAD state. 



and negative, —Be ~ 10 2 — 10 3 (in units of c 2 ). T he torus is 
described in detail in Penna, Kulkarni & Naravan (l2012l) . 

The initial magnetic field is purely poloidal. The magnetic 
field in the case of the ADAF/S ANE simulation is broken into eight 
poloidal loops of alternating polarity (Fig. [2). Each loop carries the 
same amount of magnetic flux, so the BH is unable to acquire a 
large net flux over the course of the simulation. The normalization 
of the magnetic field is adjusted such that the gas-to-magnetic pres- 
sure ratio, P, in the equatorial plane has a minimum value ^100 
for each of the eight loops. Instead of using multiple poloidal loops, 
another way of setting up an ADA F/SANE simulation is to u se a 
toroidal initial field (e.g., Model A in Igumenshche v et al .120031 and 
Model AO.OBtNIO in lMcKinnev et alJl2012l) . 

The initial magnetic field of the ADAF/MAD simulation 
forms a single poloidal loop centered at r = 300 (Fig. [3]). The 
gas accreted by the BH in this simulation has the same orientation 
of the poloidal magnetic field throughout the run, so the net flux 
around the BH increases rapidly and remains at a high value. The 
accretion flow is thus maintained in the MAD state. The minimum 
value of P in the initial torus is ~ 50. 

The magneti c field construction is described in detail in 
IPenna et alJOOl^ l 3 ! 



2.3 Preliminary Discussion of the Simulations 

The two panels in Fig. [4] show snapshots from the end of the 
ADAF/SANE and ADAF/MAD simulations. In each panel, the 
black and white streaks and red arrows show velocity streamlines 
in the poloidal plane at azimuthal angle = 0, and the dashed 
lines correspond to one density scale height. The main difference 
between the two simulations is that the SANE run exhibits more 
turbulence com pared to the MAD r un. 

Following IPenna etaD J2010h . we define the mass accretion 
rate M, the accreted specific energy e, and the accreted specific 
angular momentum j, at radius r and time t, as follows: 



M(r, t) 
e(r,t) 



-IS 

E{r,t) 
M(r, t) 

j(r,t) 
M{r,t) 



pu dAo 



TldAe 



= — II 

M(r,t) Je Jj> 

= L_ f f T r 

M(r,t) hh ^ 



(1) 

(2) 
(3) 



where dAe^ = ^f^gdOdcj) is an area element in the 0-<j) plane, 
p is rest mass density, is the four-velocity, and T[ and XJ are 



has rstart = 25M, r en d = 550M, and X B = 2.5. The ADAF/MAD 
3 In the notation of IPenna et all <2012l) , the ADAF/SANE magnetic field magnetic field has r s tart = 25M, r end = 810M, and X B = 25. 
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Figure 4. Left: Snapshot of the ADAF/SANE simulation at t = 200, 000. Black and white streaks as well as red arrows represent flow streamlines. Note the 
turbulent eddies. The blue dashed lines indicate the density scale height. Right: Snapshot of the ADAF/MAD simulation at t = 100, 000M. There is much 
less turbulence. 



components of the stress-energy tensor describing the radial flux of 
energy and angular momentum, respectively: 

T[ = (p + Tu + b 2 )u r u t -b r bt, (4) 
T; = {p + Tu + b 2 )uu^-b r b^. (5) 

The quantity u is the internal energy of the gas, Y is its adiabatic 
index which is set to 5/3 in both simulations, and b^ is a four- vector 
which describes the fluid frame magnetic field (see iGammie et al.1 
120031 for details). In equations (D-CD, the integrals are over the 
entire sphere (0 = to tt, (j) = to 2ir), and the signs are chosen 
such that M, E, J are positive when the corresponding fluxes are 
pointed inward. More useful than e is the quantity (1 — e), which 
is the "binding energy" of the accreting gas relative to infinity. 

In addition, we define 0bh to be the normalized and averaged 
magnetic flux threading e ach hemisphere of the BH horizon (see 
iTchekhovskov et alJ201lh . 

BH (*) = — L= f f \B r {r^t)\dA e ^ (6) 
2v M Jo Jcf> 

where B r is the radial component of the magnetic field and m 
is the radius of the horizon. The integral is again over the whole 
sphere, and the factor of 1 /2 is to convert the result to one hemi- 
sphere. An accretion flow trans itions to the MAD state once 0bh 
crosses a critical value ~ 50 (ITchekhovskov et all 1201 lL l2012h . 
Thus, by monitoring this quantity, we can evaluate whether a par- 
ticular simulation is in the SANE or MAD state. 

Figure|5]shows the time evolution of M, j, (1 — e) and </>bh as 
a function of time for the ADAF/SANE and ADAF/MAD simula- 
tions. The first three quantities are measured at r = 100 while the 
fourth is (by definition) evaluated at the horizon r — m. We see 
that the magnetization parameter 0b h behaves very differently in 

4 The reason for choosing r = 10 rather than r = rn is to avoid small 
deviations that sometimes arise near the horizon because of the activation 
of floors in the HARM code. Since r = 10 is well inside the inflow equi- 
librium zone at all times of interest, it is a safe choice. 
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Figure 5. Variations of M, j and (1 — e) at r = 10, and </>bh at r = tr, 
as a function of time. Solid lines correspond to the ADAF/SANE simula- 
tion and dotted lines to the ADAF/MAD simulation. Note the very different 
behaviors of the two simulations. The decrease of M with increasing time 
is explained in Fig.[6]and the text. 



the two simulations. In the ADAF/SANE simulation, 0bh remains 
small, except for one spike at time t ~ 140, 000. In contrast, in 
the ADAF/MAD simulation, the magnetization quickly rises to a 
value ~ 50 a nd remains at this high value for the rest of the run. As 
explained in ITchekhovskov et all &201 iK the plateau in 0bh cor- 
responds to the MAD state where the BH has accepted as much 
magnetic flux as it can hold for the given mass accretion rate. Any 
additional flux brought in by the accreting gas remains outside the 
horizon, where it "arrests" the accretion flow. 

Corresponding to the dramatic difference in 0bh in the two 
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Figure 6. Top Left: Shows the variation of the mean mass accretion rate 
M(r) vs r in the ADAF/SANE simulation for the six independent time 
chunks S1-S6. The colour code is as follows: SI (blue), S2 (green), S3 
(red), S4 (cyan), S5 (magenta), S6 (black). The flat region of each curve 
identifies the range of r over which the accreting gas is in inflow equilib- 
rium. This range increases monotonically with time, as one expects. Top 
Right: Similar plot for the ADAF/MAD simulation for the five time chunks 
M1-M5. Colour code: Ml (blue), M2 (green), M3 (red), M4 (cyan), M5 
(magenta). Bottom Left: An explanation for why the mass accretion rate 
shown in Fig. [5] declines secularly with time in the ADAF/SANE simula- 
tion. In each time chunk, the surface density E has to match smoothly to 
the E profile of the initial torus (dotted curve). Therefore, the decrease in 
M is purely a consequence of the initial conditions. Bottom Right: Similar 
plot for the ADAF/MAD simulation. 



simulations, there are related differences in both the binding energy 
flux (1 — e) and the specific angular momentum flux j. The quantity 
(1 — e) is about two to three times larger in the MAD simulation, 
which indicates that the MAD system has more energy flowing out 
to infinity compared to the SANE simulation. Coincident with the 
spike in B h in the ADAF/SANE simulation at t ~ 140, 000, there 
is a corresponding spike in (1 — e). During this period, the SANE 
simulation seems to have made a brief detour close to the MAD 
limit. 

The specific angular momentum flux j is about an order of 
magnitude less in the MAD simulation compared to the SANE 
simulation. Once the gas has attained the MAD state, it transfers 
very little angular momentum to the BH. Instead, angular momen- 
tum is transported out, largely through the magnetic field. This im- 
plies that an ADA F/MAD accretion flow will caus e little spin-up of 
the BH . Indeed, as Tche khovskov et all < l2Q12h and lMcKinnev et all 
d2012h have shown, if the BH has virtually any non-zero value of 
a*, an ADAF/MAD flow will cause spin-down rather than spin-up. 

Before discussing the behavior of M in Fig. \5\ we first de- 
scribe the method we use in the rest of the paper to analyze the 
time evolution of quantities. We divide the data from each simu- 
lation into a number of "time chunks" which are logarithmically 
spaced in time. In the case of the ADAF/SANE simulation we have 
six time chunks, S1-S6, with each successive chunk being twice 
as long as the previous one (Table [TJ- This logarithmic spacing is 



well-suited for the issues discussed in this paper since most of the 
quantities we are interested in show power-law behavior as a func- 
tion of both time and radius. In the case of the shorter ADAF/MAD 
simulation we divide the data into five time chunks, M1-M5 (Table 
[2]). Note that there is no overlap between chunks, and hence each 
chunk provides independent information. 

Returning to Fig. \5\ we see that M shows a large decrease 
with time in both simulations. Fig. [6| explains the reason for this. 
Since the accreting gas originates in the initial gas torus shown 
in Figs. [2] and [3] the mass distribution in the flow has to match 
smoothly to this mass reservoir. With increasing time, the radius 
range over which the flow achieves steady state increases (as dis- 
cussed in greater detail in the following sections). At the bound- 
ary of the steady state region, quantities like the surface density, 
E = (l/27r) f / pdAecf) (shown in Fig. [6]), have to match the cor- 
responding values in the torus, and this fixes M for that epoch. 
Since the torus has a prescribed variation of E with r, we thus 
have a pre-determined variation of M with time. In hindsight, it 
might have been better to design the initial torus so as to obtain a 
ro ughly constant M wi t h tim e. An alternate approach, pioneered 
by llgumenshchev et al.1 J2003), is to inject mass steadily at some 
outer radius rather than to start with a fixed total mass in a torus. 



2.4 Resolving the MRI 

Following Hawlev. Guan & Krolikl bOllh . we determine how well 
the MRI is resolved in our simulations by computing the parameters 



2tt \b e 



2tt \b*\ 



(7) 



Here, the grid cell sizes, dx d , dx^, and the magnetic field com- 
ponents, b e , b^, are evaluated in the orthonormal fluid frame. The 
fluid's angular velocity is Q. The parameter Qq is defined such that 
it becomes Amri / dz in the limit of a vertical field, where Amri is 
the wavelength of t he fas test growing mode of the linear MRI. 

lHawlevetal.l([20T 3) considered a number of diagnostics, prin- 
cipally I and dimensionless viscosity parameter a, but also 
Bz/B^ and plasma (3 = P gas /P m a S , as a function of numeri- 
cal resolution. They studied both local shearing boxes and global 
Newtonian discs and concluded that simulations with Qq > 10 and 

> 20 are sufficiently well resolved to give quantitatively con- 
verged results. They also state that simulations with smaller values 
of Q^, but correspondingly larger values of Qq, are equally good. 
Thus, we write their criterion for convergence as Q§Q^ > 200. In 

addition, they recommend that the ratio dx^ /dx r near the disc mid- 
plane should be no lar ger than 4. 



In related work, 



Sorathia et al 1 (E012h simulated global (but 



unstratified) Newtonian discs using a wide range of resolutions and 
showed that the magnetic tilt angle, which is related to the ratio 
Bl j B^ mentioned above, is a good diagnostic for evaluating con- 
vergence. On the basis of this diagnostic, they suggest that a ratio 
dx^ /dx r < 2 is sufficient for convergence, but a ratio of 4 tends to 
be somewhat under-resolved (see their Fig. 11c). Thus, their crite- 
rion is stricter than the one proposed bv lHawlev et al.ld201ll) . 

Our simulations have Qq ~ 10 — 20 throughout the ini- 
tial magnetic loops. The initial is zero because the loops are 
poloidal. For the ADAF/SANE run, the fluid inside r = 100 and 
within one density scale height of the midplane has Qq and Q^ 
between 10 — 20, i.e., Q§Q^ ~ 200, which is sufficient accord- 
ing to lHawlev et al.1 J201ll) . Our numerical grid has dx^ /dx r « 3 
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Table 1. Time Chunks in the ADAF/SANE Simulation 



Chunk Time Range (M) t chu nk/M r strict /M r loose /M 



SI 


3000-6000 


3000 


19 


23 


S2 


6000-12000 


6000 


25 


43 


S3 


12000-25000 


13000 


29 


45 


S4 


25000-50000 


25000 


43 


62 


S5 


50000-100000 


50000 


66 


92 


S6 


100000-200000 


100000 


86 


113 



Table 2. Time Chunks in the ADAF/MAD Simulation 



Chunk 


Time Range (M) 


^chunkA^ 


^strict /M 


Rloose/M 


Ml 


3000-6000 


3000 


35 


52 


M2 


6000-12000 


6000 


37 


65 


M3 


12000-25000 


13000 


69 


90 


M4 


25000-50000 


25000 


109 


128 


M5 


50000-100000 


50000 


170 


207 



at the mid-plane, which is safe according to lHawlev etal.1 (1201 ll) 
and borderline according to Is orathia et al .1 d20 1 2l) . Overall, we con- 
clude that our ADAF/SANE simulation is adequately resolved. Our 
ADAF/MAD simulation has Q$ > 100 and ~ 50, so this sim- 
ulation is very well-resolved. 

Exploring the issue of convergence further, we note that the 
gri d used in the present study is very similar to the one employed 
by iTchekhovskov et all d201 lh for simulating their MAD models. 
These authors tested convergence by increasing the number of cells 
in the cj) direction by a factor of 2, i.e., using 128 cells over the range 
(j) — 0— 2tt instead of the fiducial 64 cells. The results they obtained 
with this increased resolution agreed with those from their fiducial 
runs, indicating that 64 cells over 2ty in (j> (or 32 cells over a wedge 
of angle tt) are sufficient for convergence. Thus we are confident 
that our ADAF/MAD run ha s sufficient resolution. 

iMcKinnev etaD 120121) describe a large number of simula- 
tions, of which one sequence of models, A*BtN10, was initialized 
with a purely toroidal field. These models, which evolve into con- 
figurations similar to our ADAF/SANE simulation, used a resolu- 
tion of N r = 128, N e = 64, = 128, which is slightly dif- 
ferent from, but generally similar to , our resolution, iV r — 256, 
N e = 128, = 64. In addition. IMcKinnev et alJ l2012h con- 
sidered one high-resolution toroidal-field model, A0.94BtN10HR, 
with N r = 256, N e = 128, = 256. Looking at the detailed re- 
sults, it is not obvious that their high-resolution model is distinctly 
superior to their standard lower-resolution models. 

Based on all of the above, we believe the two simulations de- 
scribed in this paper are adequately resolved. 



3 ANALYSIS AND RESULTS 

3.1 Criteria for Convergence and Steady State 

Figure [7] shows time- averaged, 0-averaged, ^-symmetrized results 
for the final four time chunks, S3, S4, S5, S6, of the ADAF/SANE 
simulation. The strong averaging of the simulation data eliminates 
most of the turbulent fluctuations that were evident in Fig. |U and 
enables us to focus on mean properties of the flow. The accretion 
flow is geometrically thick, as expected, and the gas velocity is pre- 
dominantly inward within one scale-height of the mid-plane. At 



higher latitudes, many velocity arrows point away from the BH, in- 
dicating that there is mass outflow. At yet higher latitudes, as we 
approach the poles, the gas appears again to flow in towards the 
BH. It is therefore not obvious how much gas actually flows out to 
infinity. We discuss this question in detail in the next subsection. 

Figure [8] shows an equivalent plot for the ADAF/MAD simu- 
lation, corresponding to the final four time chunks, M2, M3, M4, 
M5. Comparing Figs.|7]and[8] the flow streamlines in the MAD run 
show more well-organized outflow behavior. There are also out- 
flowing streamlines along the axis, suggesting some kind of po- 
lar jet. However, very little energy, and practically no mass, flows 
along this jet. Therefore, for all practical purposes, the simulation 
does not have a jet. 

A critical issue for analyzing simulation data is knowing 
which regions of the solution have had sufficient time to settle down 
to a state of "inflow equilibrium", and which regions are still in the 
process of getting there. One way to do this is by looking at plots 
such as Fig. [6l and estimating "by eye" the region of steady state. 
However, a more objective criterion is preferab le, so we follow th e 
prescription for inflow equilibrium described in Pe nna et al ] {2013). 
For each time chunk, we compute the time- averaged radial velocity 
profile v r (r) of the gas within one scale-height of the mid-plane 
(the restriction to one scale-height is to enable us to focus on the 
accretion flow rather than any mass outflow or jet). From this, we 
estimate the viscous time as a function of radius r in the standard 
way: 

* visc(r) 55 KM' (8) 

We then define two criteria, one "strict" and one "loose", to es- 
timate the radius range over which the flow has achieved inflow 
equilibrium: 

tvisc (/strict) = ^chunk/2 = ttot/4, (9) 
^visc (^loose) — ^chunk — ^tot/2. (10) 

Here, tchunk is the time duration of the chunk under consideration, 
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Figure 7. Average flow properties of the ADAF/SANE simulation during chunks S3 (top left), S4 (top right), S5 (bottom left) and S6 (bottom right). In each 
panel, the flow has been averaged over the duration of the chunk £ c hunk (Table [T}, over azimuthal angle 4>, and symmetrized around the mid-plane. Colour 
indicates log p, arrows indicate direction (but not magnitude) of the mean velocity, and slanting dashed lines indicate the local density scale height. The two 
circular solid lines correspond to the steady state radius limits r str ict (thick line) and ri oose (thin line), computed using the mean radial velocity within one 
scale height of the mid-plane (see text and Table [T]f or details). 



and ttot is the total run time from the beginning of the simulation 
up to the end of the current chunl|fl 

The philosophy behind the above criteria is that we expect the 
flow to reach inflow equilibrium on a time scale of order the vis- 
cous time. Further, it takes a few viscous times to average out fluc- 
tuations. The strict criterion has ttot = 2 1 chunk = 4t v i sc , which is 
a fairly safe and conservative choice, while the loose criterion takes 
a more o ptimistic view of ho w soon inflow equilibrium is achieved. 
Note that lPenna et al.l d201Qh defined inflow equilibrium by the con- 
dition ttot = 2t v isc, which is the same as our present loose crite- 
rion. The values of t chunk, Strict and noose f° r me various time 
chunks are listed in Tables Q]and|2] and r str ict and noose are shown 
as circular solid lines in Figs.|7]and[8] It will be noticed that the ob- 
jectively determined r strict and noose are compatible with values 
one might deduce by visual inspection of Fig.[6] 

In Figs. [7] and [8] the time-averaged velocity streamlines are 
well-behaved within the respective inflow equilibrium regions of 



5 Note that the chunks are so defined that the duration of each chunk is half 
the total run time of the simulation up to that point (Tables [T] [3 



the four panels. Note also that the steady state zone is much more 
extended in the MAD simulation compared to the SANE simula- 
tion. For instance, MAD chunk M5, which has run only half as 
long as SANE chunk S6 is converged out to a substantially larger 
radius (compare the values of r s trict, noose in Tables [T] and [2]). The 
reason is the larger radial velocity of the gas in the MAD simulation 
(compare Figs. QT] and [12]). 

When the accretion flow has reached inflow equilibrium, we 
expect 6- and ^-integrated fluxes of conserved quantities, as de- 
fined in equations tH)-©, to be independent of radius. Recall that 
there is no radiative cooling, hence there ought to be strict conser- 
vation of not only mass, but also energy and angular momentum. 
As time proceeds, the range of r over which these fluxes are con- 
stant will increase, and should track r strict or noose (depending on 
the degree of constancy one requires). 

Figure [9] shows the fluxes of specific angular momentum j 
and specific binding energy (1 — e) for the six time chunks in 
the ADAF/SANE simulation. The range of radius over which these 
fluxes are in inflow equilibrium increases from time chunk S 1 to S6, 
i.e., with increasing time, as expected. The solid line segments in 
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Figure 8. Similar to Fig. [7] but for time chunks M2 (top left), M3 (top right), M4 (bottom left), M5 (bottom right) of the ADAF/MAD simulation. Note that 
in chunk M5 (lower right) r str ict and n oose both lie outside the plotted area (see the numerical values given in Table|2j. 



the plot correspond to the strict criterion r < r str ict, and the dotted 
lines correspond to the loose criterion r ^ ri OOS e- This convention 
is adopted in all later plots. 

Figure [9] highlights the difference in convergence properties 
between the two criteria. Although the strict criterion is not per- 
fect, the fluxes do remain nearly constant over the radius ranges 
defined by this criterion. The loose criterion, however, shows unac- 
ceptably large deviations from flux constancy. Hereafter, we quote 
quantitative results only for regions satisfying the strict criterion 
(the inner solid circles in Fig. |7J, though we plot results for bott0. 
Interestingly, the angular momentum flux shows larger deviations 
from constancy than either the binding energy flux (1 — e) or the 
mass accretion rate (shown in Figs.|6]and[T3j- We are not sure why 
this is the case. 

Figure [^indicates that there is a slow secular decrease in the 
converged values of both j and (1 — e) with time; the values for 
chunk S6 are smaller than those for S5, and so on. This is similar 
to, though not as extreme as, the declining trend in M already seen 



6 Obviously, more accurate results could be obtained by using an even 
stricter criterion, e.g., t v - lsc ^ ^chunk/4. However, this would reduce the 
range of r so much that we would not have sufficient dynamic range to 
obtain any useful results. 



in Fig. [6] We suspect that, in the case of j and (1 — e), the reason 
for the decline is that the SANE simulation is slowly approaching 
the MAD limit (despite our best efforts to avoid it). 

Figure [10] shows equivalent results for the ADAF/MAD simu- 
lation. Here, j and (1 — e) are less well-behaved than in the SANE 
simulation. In fact, it appears that even r s trict may overestimate the 
actual radius out to which inflow equilibrium has been achieved. 
The binding energy flux (1 — e) is a few times larger for the MAD 
simulation compared to the SANE simulation. This implies that 
the MAD accretion flow returns mechanical and magnetic energy 
to infinity more efficiently compared to the SANE simulation. In 
essence, the outflowing gas carries more energy per unit mass. The 
angular momentum flux j is substantially smaller in the MAD sim- 
ulation compared to the SANE run. Indeed j appears secularly to 
approach zero with increasing time, as seen also in the highly sub- 
Keplerian values of (compare Figs. QT] and [12]). In fact, it seems 
that BH spinup via an ADAF/MAD accretion flow is highly ineffi- 
cient. This agrees with the results reported in iTchekhovskov et all 
(120121) and lMcKinnev et alJ 120121) . 

Figure QT] shows the radial velocity |tv(r)|, the specific an- 
gular momentum i^(r) of the gas within one scale height, and the 
normalized scale height h/r. There is good internal consistency be- 
tween the profiles from successive time chunks. This is especially 
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Figure 9. The black dotted line at the top labeled jisco corresponds to 
the angular momentum of a Keplerian orbit at the radius of the IS CO. This 
represents the specific angu lar momentum flowin g into the BH in the case 
of a standard thin disc (Novikov & Thorne 1973). The cluster of lines just 
below the dotted line shows the run of specific angular momentum flux with 
radius j(r) corresponding to chunks SI (blue), S2 (green), S3 (red), S4 
(cyan), S5 (magenta) and S6 (black) for the ADAF/SANE simulation. All 
of these curves lie below the NT curve, indicating that the ADAF flow is 
sub-Keplerian, as predicted by theory. Each of the curves has a flat segment 
where the time-averaged flow shows excellent steady state convergence and 
a region at larger radii where j deviates from steady state. The bottom set of 
lines (same colour coding) shows the specific binding energy flux (1 — e) 
for the same time chunks. For both sets of lines, the solid and dotted line 
segments correspond to r ^ r str ict and r < n oose , respectively (see text 
andTablesCQEJ. 




Figure 10. Similar to Fig.|9] but for the ADAF/MAD simulation. The colour 
coding is: chunk Ml (blue), M2 (green), M3 (red), M4 (cyan), M5 (ma- 
genta). 



true when we focus only on the regions that satisfy the strict crite- 
rion for inflow equilibrium (the solid line segments). Specifically, 
apart from a tendency for h/r to increase slightly with time, the 
profiles of various quantities in successive time chunks line up well 
with one another, showing that we have a well-behaved accretion 
flow. We view the good agreement as a sign of convergence in our 
results. 

At r — 100, we have \v r \ ~ 0.002, which is far smaller than 
the local free-fall velocity w 0.14. This is to be expected. The 




Figure 11. Top Left: Shows the density- weighted mean radial velocity of 
the gas in the ADAF/SANE simulation within one scale height of the mid- 
plane during time chunks S1-S6. The colour code and line types are the 
same as in Fig.[9] Top Right: A similar plot for the density- weighted specific 
angular momentum of the accreting gas. The black dotted line shows the 
Keplerian profile of angul ar momentum for a standard thin accretion disc 
iNovikov & Thorndll973h . Bottom Left: Plot of the density scale height 
h/r for the six time chunks. Bottom Right: Plot of the mid-plane values 
of /x, which represents the normalized flux of the Bernoulli parameter (see 
eq.[T3t. The fact that is negative indicates that the mid-plane gas is bound 
to the BH. 




Figure 12. Similar to Fig. [TT] but for the ADAF/MAD simulation. Colour 
coding is as in Fig.fTo] 
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radial velocity in a viscous flow is ^ a(h/r) 2 v^, where a is the 
dimensionless viscosity parameter and (h/r) is the dimensionless 
geometrical thickness of the disc. The simulated system has h/r ~ 
0.4 and a ~ 0.05 (near r ~ 100), and this explains the observed 
velocity. 

The specific angular momentum of the accreting gas is sub- 
Keplerian (as predicted by simple ADAF models). Interestingly, u<f, 
continues to decline with decreasing radius even in the plunging 
region, i.e., inside the innermost stable circular orbit, nsco = 6. 
It appears that the dynamics of an ADAF are not strongly modified 
when the gas crosses the ISCO. This is in contrast to geometrically 
thin discs, where the angular mom entum becomes nearly constant 
once the gas flows inside the ISCO ( Shaf ee et al.l2008l ; |Penna et all 
[203). 

The fourth panel in Fig. QT] shows the normalized Bernoulli- 
flux parameter p (defined below in eq. [13]) of the mid-plane gas. 
Recall that the initial gas in the torus had Bernoulli in the range 
10 -2 — 10 -3 . The mid-plane gas in the accretion flow has a more 
negative value of p, which means it is more tightly bound to the 
BH compared to the initial gas. The profiles from the different time 
chunks agree reasonably well with one another, but not perfectly. 
This is perhaps to be expected since p is computed as the differ- 
ence of two quantities of order unity. Note that the outflowing gas 
we consider in the next subsection has a positive p. That gas has 
acquired extra energy in the pr ocess of accretion, an d it is the extra 
energy that drives the outflow dNarav an & Yi 19941) . 

Figure [12] shows the corresponding results for the MAD simu- 
lation. The radial velocity is substantially larger compared to the 
SANE simulation. Indeed, this is the reason for the larger zone 
of inflow equilibrium in this simulation. Both disc thickness h/r 
and Bernoulli Be show more fluctuations between successive time 
chunks. This is part of a pattern — fluctuations of all quantities are 
generally larger in the MAD simulation. The MAD flow is slightly 
thicker than the SANE flow, h/r ~ 0.5 compared to ~ 0.4, but it 
has roughly the same (negative) value of Be at the mid-plane. 



3.2 Mass Loss in an Outflow 

The main motivation behind the present study is to evaluate the 
amount of mass loss experienced by an ADAF through winds and 
outflows. Figures [7] and [8] show that mass does flow out in both 
the SANE and MAD simulations. However, just because a given 
parcel of gas moves away from the BH does not necessarily mean 
that it escapes to infinity. The gas might just move out for a cer- 
tain distance, turn round and merge with the inflowing gas. We 
need a physical criterion other than mere outward motion to deter- 
mine whether or not mass is lost. Before proceeding further we note 
that there is no sign of a rela tivistic polar jet in our s imulations, in 
agreement with the results of iMcKinnev et al.l d2Q12h for their runs 
with non-spinning BHs. This is perhaps not surprising since there 
is growing evidence that relativistic jets are powered by BH spin 
dTchekhovskov et alJEoill : iNaravan & McClintockll2012h . In any 
case, the discussion below is concerned with non-relativistic mass 
outflows, not jets. 

We work with gas properties averaged over the duration of a 
time chunk £ c hunk and azimuthal angle 0, and symmetrized around 
the mid-plane. We do this not only for quantities like density and 
velocity, but for all other quantities mentioned below, e.g., put, 
uut, b 2 u t , etc. As Figs. [7] and [8] show, such averaging eliminates 
all turbulent fluctuations inside the region of inflow equilibrium, 
allowing us to focus on the mean properties of the flow. This is 
important when trying to evaluate the magnitude of outflows. 



We have considered three criteria for deciding whether a gas 
streamline escapes to infinity. The first two criteria involve vari- 
ants of the Be rnoulli parame t er of t he gas. This was the parameter 
considered bv lNaravan&Yil(ll994l) in their original work in which 
they identified mass loss as being potentially important in ADAFs. 
In Newtonian hydrodynamics, Be is the sum of kinetic energy, po- 
tential energy and enthalpy. At large distance from the BH, the po- 
tential energy vanishes. Since the other two terms are positive, gas 
at infinity must have Be ^ 0. Furthermore, in steady state and in 
the absence of viscosity, Be is conserved along streamlines. Hence 
any parcel of gas that flows out with a positive value of Be can po- 
te ntially reach infinity. T his was the crux of the argument proposed 
bv lNaravan& Yil(ll994l) . 

In our case, we have an MHD flow in a general relativis- 
tic space-time. He re, the Bernoulli parameter may be written as 
dPenna et al.ll2012l) 

(pu t )+r(uu t ) + (b 2 u t ) 



Be 



(p) 



(11) 



where ( • • • ) indicates an average over time and azimuth. We sub- 
tract unity to eliminate the rest mass energy of the gas. Far from the 
BH, the expression in dTTb reduces to the Newtonian quantity — ki- 
netic energy plus gas enthalpy plus magnetic enthalpy — which has 
to be positive. Therefore, gas in a given poloidal cell of the simula- 
tion is likely to escape to infinity if the time-averaged properties in 
that cell satisfy the following two conditions: (1) the mean velocity 
has an outward radial component, i.e., (v r ) > 0, and (2) the gas 
has Be ^ 0. This is the first of three criteria we have considered. 

Because magnetic stress is anisotropic, the contribution of the 
magnetic field to the Bernoulli is not well-defined. Therefore, some 
authors (e.g., Tchekhovskoy etal.1 l201ll : lMcKinnevetai]l2012h 
ignore the magnetic term and consider the following modified 
Bernoulli parameter, 

{pu t ) + Y(uu t ) 



Be' 



(p) 



(12) 



This is arguably a more robust quantity, though it underestimates 
the Bernoulli. The second criterion we have considered for identi- 
fying outflowing gas is that it should satisfy (1) (v r ) > and (2) 

Be' ^ 0. 

Our third criterion involves a normalized energy out- 
flow rate, similar to the ratio p of energy flux to rest mass 
flux discussed in theories of magnetized relativistic jets (e.g., 
iTchekhovskov. Naravan & McKinnevll2010l) . For our general rel- 
ativistic MHD flow, we define p to be 

(If) , 



P : 



(puP) 



(13) 



where the index p refers to "poloidal", and we subtract unity to 
eliminate the contribution due to rest mass. Note that (Tf) / (pu p ) 
is just a local version of E/M in equation (0). Thus, p measures 
the flux of the Bernoulli (normalized by mass flux) and is the most 
natural quantity for our analysis. In particular, it includes the contri- 
bution of the magnetic shear stress (terms like Vb^ in eq.[5]), which 
is not included in the definitions of Be and Be' above. As before, 
we consider a parcel of gas to escape to infinity from a given radius 
r if (1) its average velocity at r is pointed outward, and (2) p ^ 0. 
For a steady axisymmetric ideal MHD flow, p is conserved along 
an outflowing streamline. Hence this /i-based criterion is arguably 
the most physically well-motivated of the three criteria, and the one 
closest in spirit to the original work of lNaravan & Yild 19941) . 

Using each of the three criteria described above, we have com- 
puted the mass outflow rate M ou t (r) as a function of r for each of 
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Figure 13. The horizontal lines near the top of the plot show the net mass 
inflow rate M(r) for the six time chunks S1-S6 of the ADAF/SANE simu- 
lation, normalized by the net mass accretion rate on to the BH, Mbh- The 
colours and line types are as in Fig. [9] The vertical lines near the bottom 
show the variation of the mass outflow rate M ou t(r) according to the p 
criterion (the results are similar to those obtained with the Be or Be' cri- 
teria), again normalized by Mbh - There is poor convergence in the results 
for the outflow, since no two successive time chunks are consistent with 
one another. The deviations are systematic — in the last three time chunks 
(S4:cyan, S5:magenta, S6:black), each successive time chunk gives a lower 
M ut at a given r compared to the previous chunk. Hence, the mass outflow 
rates shown here should be interpreted as upper limits. 
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Figure 14. Mass outflow rate in the ADAF/MAD simulation based on the 
p criterion. The colours and line types are as in Fig. [10] The last three 
chunks (M3:red, M4:cyan, M5:magenta) show large and systematic devia- 
tions, suggesting that (as in the case of the ADAF/SANE simulation) we do 
not have good convergence and the computed mass outflow estimates are 
upper limits. 



the time chunks in the ADAF/SANE and ADAF/MAD simulations. 
The results from the three criteria agree well with one another. We 
show plots corresponding to only the p criterion. 

Figure [13] shows for the ADAF/SANE simulation the mass 
outflow rate M ou t (r) and the net mass inflow rate M(r), both nor- 
malized by the net mass accretion rate on to the BH, Mbh. The 
results for the mass inflow rate M(r) are identical to those shown 
in the top left panel of Fig. [6] except that the normalization by Mbh 



shifts the curves vertically and causes them to lie on top of one an- 
other. 

Surprisingly, the results for M ou t show very poor conver- 
gence. Specifically, the M ou t profiles corresponding to different 
time chunks deviate substantially from one another. Moreover, the 
deviations are systematic. In each time chunk, the outflow appears 
to pick up just around the limiting radius for inflow equilibrium. 
Since the latter moves out for later chunks, the entire M ou t profile 
also moves out. Apparently, at each time, the current estimate of the 
mass outflow rate at a given radius is an overestimate compared to 
the rate we would estimate at a later time (compare in particular the 
last three time chunks shown in cyan, magenta and black). Because 
of this, the outflow rate estimate even from the last time chunk S6 
(black curve) must be viewed only as an upper limit. Moreover, 
even this estimate corresponds to a mass loss rate at r ~ 100 no 
more than the net inflow rate Mbh into the BH. Given that it is an 
upper limit, we can state with some confidence that mass outflow 
is unimportant for m < r < 100. 

It is useful to co mpare our results with those obtained by 
iMcKinnev etaD l2012h for their model AO.OBtNlO. This model 
was initialized with a toroidal field and is an excellent example of 
an ADAF/SANE system. In Table 4 of their paper, the authors pro- 
vide various estimates of the mass outflow rate measured at a char- 
acteristic radius r G = 50. Their quantity M mw ,o is most relevant 
since it focuses on unbound gas, defined as Be' > OB The normal- 
ized mass outflow rate, M m w,o/MH, that IMcKinnev et aD 120121 ) 
find at r = 50 in model AO.OBtNlO is essentially zero, in good 
agreement with our result, M ou t/MBH = 0.07 at r — 50 in chunk 
S6; at r = r s trict = 86, our outflow rate is M out / Mbh = 0.6. 
It should be noted that M m w,o includes additional constraints, viz., 
that the escaping gas should have b 2 / p < 1 and gas to magnetic 
pressure ratio f3 < 2. Our mass outflow criteria do not include 
these constraints. When we include them, we find that our mass 
outflow rate is zero at r = 50 and 86. Apart fr om these details, 
both t he present work and model AO.OBtNlO in Mc Kinnev et al.l 
(120121) agree on the following key result: out to radii ~ 50 — 100, 
ADAF/SANE systems have negligible mass outflow. 

Figure Q3] shows mass outflow estimates obtained via the p 
criterion for the ADAF/MAD simulation. As in the case of the 
ADAF/SANE simulation, the convergence behavior is poor. In par- 
ticular, the results from chunks M3 (red), M4 (cyan) and M5 (ma- 
genta) do not agree well with one another. Thus, once again, we 
believe the mass outflow rates we estimate from this simulation 
should be viewed as upper limits. 

Despite the unsatisfactory convergence, if we take the results 
at face value, we find for time chunk M5, M ou t / Mbu ~ 0.2, 0.6, 
1.1, at radii r = 50, 100, 170 (= r s trictX respect ively. Two of 
the simulations described in IMcKinnev etaO l2012h . AO.OBfNIO 
and AO. ON 100, correspond to MAD flows around non- spinning 
BHs and are good comparisons (though our simulation has run 
significantly longer). At radius r Q = 50, AO.OBfNIO has essen- 
tially zero outflow, i.e., M mW)0 /MH ~ 0, while A0.0N100 has 
Mmwo/M H « 0.4. Our estimate, M ou t/MBH ~ 0.2, agrees 
welfl 

7 The authors define a second quantity, M W)Q , which represents all out- 
flowing gas, regardless of whether the Bernoulli is positive or negative. It 
is less relevant for us since most of this gas is bound to the BH and can- 
not escape to infinity. We thank J. McKinney (private communication) for 
clarifying the definitions of M mw , and M w . . 

8 As mentioned earlier, McKinn ev et al.l l2012h require several conditions 
to be satisfied, viz., v r > 0, Be' > 0, b 2 /p < 1, /3 < 2, before they 
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Figure 15. Left: Shows five outflowing streamlines in time chunk S6 of the 
ADAF/SANE simulation. The streamlines have their footpoints at (r, 6) = 
(86,0.2), (86,0.4), (86,0.6), (86,0.8), (86,1.0). All five streamlines 
have positive values of fi at their footpoints. Right: The variation of /i along 
each of the streamlines in the left panel, using the same line types. Note that 
fi shows large deviations from constancy for the last two streamlines. 
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Figure 16. Similar to Fig. Q2] but for the ADAF/MAD simulation. The 
streamline footpoints are at (r, 9) = (170,0.2), (170,0.4), (170,0.6), 
(170, 0.8). All four streamlines have positive /i at their footpoints, and all 
show good conservation of [i. 



We have looked a little deeper into why the M ou t (r) profiles 
we obtain from our simulations show poor convergence. Figure [15] 
shows results corresponding to five streamlines in time chunk S6 
of the ADAF/SANE simulation. These streamlines have footpoints 
at r = rstrict = 86 and = 0.2, 0.4, 0.6, 0.8, 1.0 rad, respec- 
tively. All these streamlines have a positive value of \i at their foot- 
points. Since fi is supposed to be conserved along each streamline, 
all of this gas ought to escape. The right panel of Fig. Q3] shows 
the variation of \i along each streamline as the gas moves away 
from the BH. We see that \i is approximately constant and positive 
for the the three streamlines closest to the pole. However, the two 
streamlines closer to the disc show a sudden drop in the value of \i 
as one moves outward. Clearly these streamlines have not reached 
steady state, since n would then be constant. It seems likely that 
the positive value of \i for these streamlines is a transient feature. 
Unfortunately, these suspect streamlines carry the most mass. 

Figure [16] shows similar results for four outflowing stream- 
lines in the ADAF/MAD simulation. Here the conservation of \i 
along outgoing streamlines is satisfied much better. In addition, the 



include a particular gas streamline in their estimate of M mW)0 . When we 
apply the same conditions on our ADAF/MAD simulation, we estimate the 
mass outflow rate at r = 50 to be 0.06, still in good agreement with their 
outflow rates. 
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Figure 17. Analysis of convective stability of the ADAF/SANE simulation. 
Results are shown for time chunk S6 using time- and azimuth-averaged, 
symmetrized, simulation data. At each point (R, z) in the poloidal plane, 
the maximum growth rate 7 according to the two Hoiland criteria are com- 
puted. Stable regions are shown by blank areas. Unstable regions where 7 < 
Qk/30 are indicated by crosses, regions with f2x/30 ^ 7 < Q^/10 are 
indicated by open circles, and regions with 7 ^ Qr/ 10 are indicated by 
filled circles. The solid and dotted lines correspond to one and two density 
scale heights, respectively. Note that the accretion flow is stable to con- 
vection over the entire inflow region. The instability near the poles is not 
significant since the analysis is not valid there. 



value of fi is generally larger, which indicates that the outflowing 
gas carries more energy per unit rest mass. 



3.3 Convection 

A secondary goal of this study is to investigate the importance of 
convection in magnetized ADAFs. It is well-known that the entropy 
profile in an ADAF has a large negative gradient, making the flow 
highly unstable by the Schwarzschild criterion. However, an ADAF 
also has angular momentum increasing outward, which has a stabi- 
lizing effect on convection. 

For axisymmetric rotating flows, the two Hoiland criteria de- 
termine whether or not gas is convectively unstable. The same 
criteria are likely to remain approximately valid also in magne- 
tized flows, so long as the field is reasonably weak, since the 
long- wavelength con vective modes are effectively hydrodynamical 
dNaravan et al.l l2002). In addition, since convection is alocal insta- 
bility, the relativistic versions of the Hoiland criteria dSeguinl 19751) 
carry over directly to general relativity by the equivalence principle. 

We have analyzed the final time chunk S6 in the ADAF/SANE 
simulation to determine the level of convective instability in the ac- 
cretion flow. Figure[T7]shows the result. In brief, all the fluid within 
two scale heights of the mid-plane appears to be convectively sta- 
ble. The gas is certainly turbulent (see Fig.|4j - this is what enables 
it to accrete - but it is apparently not convective, at least by the Hoi- 
land criteria. Rather, the turbulence seems to be entirely the result 
of the MRL Could magnetic fields be confusing the issue? We think 
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Figure 18. Similar to Fig.QTJbut for the ADAF/MAD simulation. 



this is unlikely. An alytical studies of convection in the presence 
of magnetic fields (Balb us&Hawievll2002l : | Naravan et al 
show that magnetic fields generally act in such a way as to stabilize 
convection. That is, a fluid configuration that is convectively unsta- 
ble could be made stable by a suitable field, but not the other way 
round. Of course, the magnetic field might induce its own instabil- 
ity, e.g., MRI, but this can no longer be considered convection. We 
intend to explore this question in greater depth in the future. 

Figure[l8]shows the convection properties of the ADAF/MAD 
simulation. Based on the Hoiland criteria, it appears that the 
MAD simulation is more unstable to convection compared to the 
ADAF/SANE simulation. This is not surprising. The gas rotates 
much more slowly and hence the stabilizing effect of rotation, 
which we think is the primary reason for the lack of convection 
in the ADAF/SANE simulation, is no longer effective. We caution, 
however, that the magnetic stress is larger in the MAD simulation, 
and the Hoiland criteria do not include the effect of this stress. By 
the argument in the previous paragraph, the magnetic field might 
well be strong enough to switch off the convective instability even 
in those regions where the Hoiland criteria indicate instability. The 
accreting gas in the MAD simulation has very little turbulence, so 
it certainly does not manifest any of the usual features of turbu- 
lent convection. We suspe ct that the flow i s in a state of frustrated 
convection as proposed bv lPen et al ]J2003|). 



4 ADAF OR CDAF OR ADIOS? 

As originally defined, an ADAF is any accretion flow in which en- 
ergy advection is more important than energy loss through radia- 
tion. In this sense, the term is all-inclusive. However, sometimes 
the name ADAF is used in a more restrictive sense, where the flow 
is not only advection-dominated but also has negligible mass loss 
through a wind and is not strongly convective. If we further restrict 
ourselves to a flow that shows self- similar behavior, we have the 



classic ADAF scalings dNaravan & Yilll994l : lNaravan et al.|[l998L 

v r ~-ar- 1/2 , p~Mor\-*l\ Q ~ (5/3 - r) 1/2 r _3/2 , 

(14) 

where M is the steady mass accretion rate, a is the viscosity param- 
eter, Q is the angular velocity, and V is the adiabatic index. These 
scalings follow from basic conservation laws and the a prescription 
for viscosity. By assumption, there is no mass outflow. 

In the same spirit, the convection-dominated accre tion flow 
(CDAF, Naravan et al. 2000; Ouataert & Gruzinov 2000bl) is an ac- 
cretion flow in which the dynamics are determined by conservation 
laws plus a steady outward flux of energy carried by convection. 
This requirement gives the following CDAF scalings, 

-3/2 -1/2 -3/2 /1C x 

v r ~ —r 1 , p ~ Mr 1 , il ~ r 1 . (15) 

Once again, there is no mass outflow. 

Final ly, the advection-dominated inflow-outflow solution 
(ADIOS, Blandford & Begelman 1999) describes a system in 
which a strong wind carries away mass, angular momentum and 
energy. Nothing is conserved in this model, so there is consider- 
able freedom in the form of the solution. It is generally assumed 
that quantities behave as power-laws of radius, which motivates the 
following ADIOS scalings, 

-1/2 -3/2+s -3/2 

v r ~ —ctr , p~r ' , \l ~ r , (16) 

where s is a free index which can have a value anywhere between 
(self-similar ADAF) and 1 (maximal ADIOS). The mass outflow 
rate in this model scales as M ou t oc r s . Recently. Ib egelmar] d2012h 
has presented arguments suggesting that s ~ 1. 

All of the above models are based on a fluid description, with- 
out allowing explicitly for magnetic fields. We believe this is rea- 
sonable, at least for the ADAF/SANE simulation, where the mag- 
netic stress behaves to a good approximation like viscosity, and the 
magnetic pressure is no t very important relative to gas pressure. 
lAkizuki & Fukud (120061) have developed self- similar solutions for 
magnetized ADAFs. However, they assume a purely toroidal field 
(no shear stress) and consequently have to invoke a-viscosity. 
Moreover, their solutions are similar to the ADAF/ ADIOS solu- 
tions mentioned above so long as the magnetic pressure is modest, 
as in the ADAF/SANE simulation. This last condition may not be 
true for the ADA F/MAD simulation. How ever, even for a MAD 
flow, the model of A kizuki & Fukuel j2006h is not appropriate since 
it assumes a toroidal field, whereas the key feature of the MAD so- 
lution is a strong poloidal field. 

We have shown in |3]that the ADAF/SANE and ADAF/MAD 
simulations appear not to be convective, to the extent we can tell 
from the Hoiland criteria. We did not include the effect of the mag- 
netic field, so we cannot make any firm statements regarding con- 
vection. Nevertheless, for the present, we will assume that neither 
simulation is a full-fledged CDAF. Also, neither flow has signifi- 
cant mass outflow up to r ~ 100. We can thus say that the simula- 
tions are best described as "basic" ADAF^l over this radius range, 
though it is possible that they are just beginning to make a transi- 
tion to the ADIOS state beyond r = 100. From equations dT4l) and 
(Q2}, we see that both solutions predict \v r \ ~ or -1//2 , which can 
be checked. 

The left panel in Fig.[T9lshows the velocity profiles in the final 

9 By "basic ADAF" we simply mean an ADAF that has no convection 
and no significant outflows. Systems with convection (CDAFs) and strong 
outflows (ADIOS) are still ADAFs in the general sense of the term, but they 
are not "basic ADAFs". 
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Figure 19. Left: Radial velocity |^r(V)l °f th e § as i n ti me chunk S6 of 
the ADAF/SANE simulation (from Fig. [TT} and time chunk M5 of the 
ADAF/MAD simulation (from Fig.fl2l. The two dashed lines have slope 
equal to — 1 /2, the value expected in the self-similar regime for both a basic 
ADAF and an ADIOS. Over most of the volume, the velocity varies more 
rapidly with radius than expected for a self-similar solution. Right: Similar 
to the previous panel, but showing the quantity \v r (r)\/a(r). Note that the 
ADAF/SANE model agrees much better with the self-similar model, except 
as the gas approaches the ISCO (nsco — 6). 



time chunks, S6 and M5, of the ADAF/SANE and ADAF/MAD 
simulations. There is some indication that, at the outermost radii 
of the respective converged regions, the velocity is settling to the 
expected r -1 ^ 2 dependence. However, over most of the flow, the 
velocity varies more steeply with radius. Part of the explanation is 
that, in the self- similar regime, the radial velocity of an ADAF is 
approximately given by \v r \ ~ a{h/r) 2 v^ ~ 10 _2 ^ff. However, 
at the BH horizon the gas must have \v r \ = v& = c. The radial 
velocity thus has to transition from its self- similar value to the free- 
fall value. It takes a substantial range of r to achieve this, especially 
in the ADAF/SANE simulation. The radial velocity in the MAD 
simulation is larger, \v r \ ~ O.lfff, so this flow is able to follow the 
self-similar scaling closer to the BH. 

A second effect is also in operation, viz., the effective a of 
the accreting gas varies with r. The right panel in Fig. [T9l corrects 
for this by plotting \v r \/a, where a(r) is estimated directly from 
simulation data for gas within one density scale-height of the mid- 
plane. The ADAF/SANE simulation now shows satisfactory self- 
similar behavior over a wider range of r. Removing the a scaling 
does not improve things much for the ADAF/MAD simulation. 

All of this discussion is based on the radial velocity v r (r), 
which we feel is the natural dynamical variable to consider. 
Most previous authors have focused instead on the density pro- 
file p(r). In steady state the two quantities are simply related: 
M rsj pv r r 2 (h/r) ~ constant. The mid-plane density profiles in 
our two simulations are roughly compatible with the velocity re - 
sults shown in Fig. [19] Many authors, no tably lYuan et all J2012bh . 
find that the density follows a single power-law over a wide range 
of radius. The velocity does not show this property (Fig. [19). 

Figure [20] shows the dependence of the gas angular velocity Q, 
in our two simulations. The ADAF/SANE simulation shows excel- 
lent convergence in the sense that the Q(r) curves from different 
time chunks agree very well with one another. Moreover, the an- 
gular velocity follows the analytical r~ 3 ^ 2 scaling quite accurately. 
However, the normalization is not correct. Since r = 5/3, the self- 
similar ADAF model predicts Q ~ (see eq.ll4l. whereas we find 
distinctly non-zero rotation in our simulation. 

The likely explanation is that the s imulation behaves, no t like 
the steady state self- similar solution of iNarayan" & Yil dl994l) . but 



Figure 20. Left: Angular velocity Q(r) of the gas in time chunks S1-S6 of 
the ADAF/SANE simulation. The dashed line has a slope equal to the self- 
similar value of —3/2. Right: Similar plot corresponding to the five time 
chunks M1-M5 of the ADAF/MAD simulation. 




the steady state selt-similar solution or JNarayan &l Yi (1994), but 
rather like the similarity solution derived by lOgilviel (1 19991) . The 



Figure 21. Left: Radial velocities vs r for time chunks S1-S6 of the 
ADAF/SANE simulation. The colour code is the same as in Fig. [9] Note 
that each curve dives down suddenly at a certain radius. This is the stag- 
nation radius for that time chunk. Beyond this radius, the mean velocity is 
outward because of the viscous relaxation of the initial torus. Right: Cor- 
responding results for the ADAF/MAD simulation, with colour code as in 
Fig. [10] 



latter solution describes the evolution of an advection-dominated 
flow as a function of both r and t, starting from an initial narrow 
ring of material. With increasing time, the flow evolves in a self- 
similar fashion. Most interestingly, in Ogilvie's solution, the an- 
gular velocity does not go to zero anywhere except in the region 
r — > 0. In fact, over most of the volume, the rotation rate remains 
a substantial fraction of the Keplerian rate, exactly as in our sim- 
ulations. Since we started our simulations with an initial torus of 
material, the similarity solution is a better point of reference than 
the self-similar solution; the latter is valid only at asymptotically 
late time when the flow has reached steady state at all r. 

As a further comparison between the ADAF/SANE simulation 
and Ogilvie's (1999) similarity solution, Fig.[2T]displays again the 
radial velocity profiles for different time chunks, but now shown 
over an extended range of radius. The velocity in each profile 
dives suddenly to zero and becomes negative at a "stagnation" ra- 
dius r st ag. We see that r s t ag increases with increasing time, as ex- 
pected for the similarity solution. The analytical solution predicts 
t 2/3 , which means that r s tag should increase by a factor 
~ 10 between chunks SI and S6. The actual increase is a factor of 
20. We view this as good agreement. 

The ADAF/MAD simulation results shown in the right panels 
of Figs.[20]and[2T]are less convincing. This simulation has a strong 
magnetic field and an arrested mode of accretion which, based on 
the evidence of all the diagnostics plotted in various figures, makes 
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the flow behave more erratically. Analytically, the MAD regime is 
sufficiently different from the SANE regime that we cannot expect 
either the self- similar ADAF solution or Ogilvie's similarity solu- 
tion to be a good description. 

As already stated, there is a hint near the outer edges of the 
ADAF/SANE and ADAF/MAD simulations that ADIOS-like be- 
havior is beginning to take hold. If we had a larger range of radius 
in inflow equilibrium, it might be possible to estimate how the out- 
flow rate varies with radius and thereby determine the index s in 
the scaling M ou t oc r s . Unfortunately, t his is out of reach with our 
current simulations. lYuan et aD d2012bh estimate from their large 
dynamic range 2D hydrodynamic simulations that s ~ 0.4 — 0.5. 



5 SUMMARY AND DISCUSSION 

The main highlights of the present work are: (1) We have run our 
simulations for an unusually long time in an effort to approach a 
steady state ADAF as closely as possible over a wide range of ra- 
dius. (2) We have explored the role of the initial magnetic field 
topology. With respect to the latter, we have considered two very 
different limits: (1) an ADAF/SANE simulation (SANE = "stan- 
dard and normal evolution"), which is a good proxy for an ADAF 
model in which the magnetic field is merely an agent that causes 
angular momentum transport ("viscosity") but plays no important 
dynamical role, and (2) an ADAF/MAD simulation (MAD = "mag- 
netically arrested disc"), where the magnetic field is strong enough 
to alter substantially the dynamics of t he gas and to drive the sys- 
tem to a magnetically arrested state fe umenshc hev et al.l 12 003: 



magn etica lly 
Nara van et al.l 2003; Tchekhovskov et al. 
20121) . 



2011 



McKinnev et al. 



Our key result is that, for radii out to r ~ 100 (gravitational 
units, GM/c 2 ), there is not much mass loss to an outflow. Tur- 
bulence certainly leads to both inward and outward gas motions. 
However, when we consider the time- averaged gas flow and how 
much gas flows out with enough energy to escape from the gravita- 
tional potential of the BH, it turns out to be only a fraction of the net 
mass accretion rate Mbh into the BH. Quantitatively, at r ~ 100, 
we find Mout ~ 0.6 Mbh for both simulations. Furthermore we 
view these estimates as upper limits since the simulations reveal 
poor convergence in M ou t (see Figs. [T3l[T4l) . 

Because of the very long run times of our simulations, we are 
unable to run multiple realizations of the SANE and MAD configu- 
rations to explore variability from one realization to another. On the 
other hand, the long run time allows us to explore convergence as a 
function of time within each simulation. We do this by dividing the 
simulation data into a number of independent chunks in log t ( 32.31 
and Tables [T] [3. By comparing different time chunks and checking 
how any quantity of interest varies from one chunk to the next, we 
are able to decide how reliable the results are for that quantity. 

A second important issue is the range of r over which each 
time chunk has reached inflow equilibrium. We use two different 
criteria, a strict one (eq.|9) and a loose one (eq.[l0j, and estimate for 
a given chunk the limiting radii, r str ict and noose, corresponding to 
each of these criteria (Tables Q] [2). Many properties of the gas show 
good convergence among different time chunks when we limit our 
attention to radii r < r strict- The results are less convincing with 
the loose criterion. However, even with the strict criterion, we find 
that some questions such as the amount of mass loss in outflows 
cannot be answered with confidence. 

We initialized the ADAF/SANE simulation with a number of 
poloidal magnetic loops (Fig. [3 in an attempt to achieve an ac- 



cretion flow with very little net flux at each radius. By and large 
this simulation behaved the way we hoped it would. In particular, 
the magnetic flux at the BH horizon, measured by the parameter 
0bh, did not come close to the limiting MAD value of 50 (except 
for one brief glitch at time t ~ 140, 000, see Fig. [5]). Thus we be- 
lieve the ADAF/SANE simulation is a believable representation of 
an ADAF system. We could have avoided the MAD regime more 
effectively by st arting the simulation with a purely toroidal field, 
as in Model A of llgumenshchev et all d2003h or Model AO.OBtNIO 
of iMcKinnev et al j 120121) . This option is worth exploring in the 
future. 

The ADAF/SANE simulation shows good convergence and 
behaves as expected. The radial velocity, angular velocity, angu- 
lar momentum and disc thickness profiles as a function of r agree 
well between different time chunks (Figs.[TT][20]). At large radii, the 
radial velocity falls well below free-fall (Fig. [19). This is expected 
since accretion is mediated by "viscous" angular momentum trans- 
port which causes the velocity to be suppressed by a factor of a 
relative to free-fall; there is also a factor of (h/r) 2 which causes 
a further decrease in the velocity. Interestingly, as discussed in §|4] 
the ADAF /SANE simulat ion is better described by the similarity 
solution of lOgilviel dl999h than the original self- similar solution of 
iNaravan & Yil dl994l) . Nevertheless, the radial dependence of ve- 
locity follows the self-similar solution quite well (Fig. \\9\ right 
panel). 

The ADAF/MAD simulation shows quite different behavior 
compared to the ADAF/SANE simulation. The inflow velocity is 
substantially larger and the angular momentum and angular ve- 
locity are substantially smaller (Figs. [12] [20). The latter appears 
to be an important charact eristic of MAD flows. As discussed in 
Tchekhovskov et al the gas brings in very little angular 

momentum to the BH and therefore induces little spin-up even for a 
non- spinning BH. In the case of a spinning BH, a MAD flow actu- 
ally causes spin-down. The reduced rotation rate of the gas means 
that there is less centrifugal support. Consequently, the radial dy- 
namics are dominated by balance between gravity, gas pressure and 
magnetic stress. We find that the gas accretes at about a tenth of the 
free-fall speed, which is a factor of several larger than the velocity 
in the ADAF/SANE simulation. 

Because of the larger radial velocity, the ADAF/MAD simu- 
lation reaches inflow equilibrium over a substantially larger range 
of radius at a given time relative to the ADAF/SANE simulation 
(compare Tables [2] and [TJ- On the other hand, convergence in the 
sense of agreement between different time chunks is less convinc- 
ing. We suspect that the reason is the large-scale ordered magnetic 
field in the MAD simulation, which imposes coherent long-lived 
structure in the flow. 

In terms of the amount of mass outflow, the ADAF/SANE and 
ADAF/MAD simulations behave rather similarly. We tried three 
different criteria to determine how much gas escapes to infinity at 
a given radius: one criterion was based on the Bernoulli parame- 
ter Be (eq.QT|), a second on a different Bernoulli Be that ignores 
the magnetic contribution (eq. [12]), and a third on the normalized 
energy flux \i (eq. [13]). The results are nearly identical with all 
three criteria, which is reassuring. Unfortunately, the results show 
poor convergence with time. In particular, the radial variation of 
Mout(r) for the last few time chunks (S4-S6 and M3-M5) differ 
by much more than we would expect for a converged simulation. 
Nevertheless, taking the results at face value, we conclude that the 
mass outflow rate M ou t becomes comparable to the net inflow rate 
Mbh into the BH at a radius r ~ 120 = 60r H in the ADAF/SANE 
simulation and r ~ 160 = 80rH in the ADAF/MAD simulation. 
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These radii are fairly far from the BH. In fact, since our mass out- 
flow rates are upper limits, the critical radii where mass outflows 
begin to dominate could be substantially larger. 

Our result that outflows are weak out to r > 100 disagrees 
strongly with previous work. Many simulations of ADAFs have 
been described in the literature (see §1 for a brief review), and most 
of these studies have concluded that there are powerful mass out- 
flows at radii well below r = 100. On investigation, it appears that 
there is a significant methodological difference between our ap- 
proach and that used by previous authors. As explained in 33.21 all 
of our calculations are based on time- and azimuth- averaged quan- 
tities in which fluctuations due to turbulence have been eliminated. 
Only if the average velocity of gas in a grid cell has a positive ra- 
dial component, and furthermore if the gas has enough energy to 
escape from the system (/i > 0), do we consider the particular gas 
packet to be part of an outflow. Most other authors have focused 
on individual snapshots of their simulations and counted any gas 
that happened to be moving away from the BH as outflow. Since 
turbulence causes gas to move to and fro, a good fraction of the 
gas in any snapshot would be moving out simply as part of turbu- 
lent eddies. However, very little of this gas would actually leave the 
system since the velocity vector is likely to turn round on an eddy 
time. Moreover, much of the gas would probably have insufficient 
energy (/i < 0) to climb out of the BH potential. Indeed, several 
previous authors have noted, after presenting very large estimates 
for the mass outflow rate, that most of the gas in their "outflows" 
has a negative Bernoulli. 

The distinction between the approach taken in previous papers 
and in the present work can be appreciated by comparing Fig. |4] 
and Fig. [T] The snapshot of the ADAF/SANE simulation in the left 
panel of Fig. |4] shows turbulent eddies down to quite small radii. 
A fraction of the gas in each of these eddies is temporarily moving 
outward, but none of it is likely to escape to infinity. However, in 
the standard approach used to estimate the mass outflow rate, the 
outward-moving part of each eddy would be included as part of 
M ut • This is likely to lead to a large overestimate of the mass out- 
flow rate. In contrast, our calculations use the average flow stream- 
lines shown in Fig. [7] Consider the final time chunk S6 (lower right 
panel). Inside r ~ 30 — 40, there are no streamlines with velocity 
vectors pointed away from the BH. Therefore, when we compute 
the mass outflow rate, we obtain vanishingly small values of M out 
for radii < 30 (Fig. US. 

Because of the above major difference between our calcula- 
tions and those of pre vious authors, it is hard t o compare our results. 
The one exception is iMcKinnev et"aD l2012h . who, though basing 
their work on snapshot data, explain their calculations in sufficient 
detail to enable a comparison. Leaving aside jets, which are not rel- 
evant for the non-spinning BHs considered here, IMcKinnev et al.l 
present two distinct estimates of the mass outflow rate. One 
estimate is called M mw , and it focuses on outflowing gas with pos- 
itive Be' (it also imposes a couple of other constraints, see 33. 2t . 
This quantity is closest to our prescription for estimating the mass 
outflow. Their second outflow estimate is called M w , and it in- 
cludes essentially all outflowing gas in each snapshot, independent 
of Be. This quantity is close in spirit to mass outflow estimates in 
many other papers in the literature, and is in our view an overesti- 
mate of the actual mass loss rate because it includes gas churning 
in turbulent eddies. 

For their Model AO.OBtNIO, which is an excellent ex- 
ample of an ADAF/SA NE system around a non-spinning BH, 
IMcKinnev et alJ l2012h estimate M W /M H ~ 1.2 at r = 50 
(here Mh is the net mass accretion rate into the BH, similar to 



our Mbh), which suggests a strong outflow already at this ra- 
dius. However, they find M mw /MH to be essentially zero. In our 
ADAF/SANE simulation, at r = 50 we find M ou t/M B H = 
0.07, i.e., practically zero, in good agreement with M mw . In the 
case of their two ADAF/M AD systems around non - spinning BHs, 
AO.OBfNIO and A0.0N100, McKi nnev et aD l2012h find at r — 50 
that Mmw/ Mh = 0, 0.4, and M W /M H = 0.6, 1.1, respectively. 
Our ADAF/MAD simulation gives M ou t/M B H = 0.2, in agree- 
ment with the M mw estimates. It thu s appears that our result s are 
perfectly compatible with the work of McKin nev et aD b012h . We 
are also in agreement with lPang et al.ld201ll) « though the latter work 
is mostly concerned with the accretion of slowly-rotating gas. 

Some papers have argued for strong outflows based simply on 
the fact that the radial profile of density and/or velocity do not fol- 
low the standard ADAF scalings given in equation dT4l >. Focusing 
on the radial velocity, the simulations generally show \v r \ increas- 
ing more rapidly with decreasing radius than expected in the self- 
similar solution. Our simulations too show this effect (Fig. [T9}. It 
turns out that two separate effects, neither involving outflows, cause 
the velocity profile to be modified. 

First, because the accreting gas makes a sonic transition as it 
approaches the BH and switches to a free-fall mode inside this ra- 
dius, we have \v r \ ~ vs near the BH. However, the velocity in the 



self- similar regime is far below free-fall: 



a(h/r) 



The 



flow needs a considerable range of r to adjust from one scaling to 
the other, and we believe this is a large part of the reason why the 
velocity profiles seen in simulations look so different from the sim- 
ple power-law given in equation dl4l . Cl ear examples of this effect 
may b e seen in the global ID models of iNaravan. Kato & Honmal 
dl997h , where the non- self- similar zone extends from the inner 
boundary to a few tens of gravitational radii. 

Secondly, it is the quantity v r /a that is expected to be self- 
similar, not v r itself. Since a varies with radius in our simulations 
(especially in the ADAF/SANE simulation), this causes an addi- 
tional deviation in v r (r). As Fig. [19] shows, removing the a depen- 
dence gives a better-beha ved velocity profile th at agrees fairly well 
with the models shown in lNaravan et alJ dl997h . 

Another argument for strong outflows that is sometimes used 
in the literature is to take the gas density at the outer radius of the 
simulation, and to calculate from it the Bondi mass accretion rate 
Mb - If the actual mass accretion rate Mbh into the BH in the sim- 
ulation is much smaller than Mb , then it is claimed that the differ- 
ence is because most of the incoming gas was ejected in an outflow. 
The problem with this argument is that, for a given outer boundary 
condition on the density, theory says that the accretion rate via an 
ADAF will be smaller than Mb by a factor ~ a(h/r) 2 ~ few %. 
Thus, having Mbh <C Mb is perfectly natural for an ADAF; it 
does not imply strong outflows. Note, however, that this explana- 
tion only goes so far. If it turns out that Mbh is much smaller than 
even a(h/r) 2 Mb, then one has to look for other explanations such 
as strong outflows or convection. To our knowledge, no simulation 
to date has come close to violating this limit. 

ADAFs in nature usually extend over many decades in radius. 
The ADAF around Sgr A*, for instance, extends from the BH out 
to the Bondi radius at r > 10 5 . Supermassive BHs in other low- 
luminosity AGN similarly have ADAFs extending over 5 or more 
decades in radius. In the case of stellar-mass BHs in X-ray binaries, 
the AD AF is usually formed by evapor ation from a thin disc on the 
outside dNaravan & McClintockl 120081) . For systems in quiesence, 
where the mass accretion rate is low, the transition radius is typi- 
cally ~ 10 3 — 10 4 . In contrast, simulations of ADAFs are gener- 
ally restricted to a much smaller range of radius (but see the recent 
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work of I Yuan et al.ll2012bh . How relevant are simulation results to 
real systems? 

Our views on this question are driven by insights gained 
fr om global ID mode ls o f ADAFs such as the ones show n 
in lNaravan et al.1 dl997h and lChen. Abramowicz & Lasotal d 19971) . 
These solutions show three zones: an inner zone where the flow ad- 
justs to the free-fall boundary condition at the BH, an outer zone 
where it adjusts to whatever outer boundary condition is present in 
the system (Bondi or disc evaporation), and a middle zone where 
the flow is more or less self-similar. If a simulation covers a large 
enough radius range to capture some piece of the middle zone, then 
it would be straightforward to stretch out the self- similar regime to 
any radius range we require. We suspect that the two simulations 
presented in this paper may have just managed to develop a piece 
of the middle zone, but we do not have any proof of this. In any 
case, we believe that only by obtaining inflow equilibrium over a 
sufficiently large range of radius can we hope to use simulations to 
make useful statements about real flows. 

It should be noted that the properties of the self-similar mid- 
dle zone are fairly insensitive to parameters. There is an obvious 
dependence on a (see eq. [T4]) and a modest dependence on rF°l 
but virtually nothing else matters. In other words, provided ADAF 
conditions are satisfied, the accretion flow will head towards the 
particular disc thickness h/r and Bernoulli Be(r) it wants in the 
middle zone, regardless of the precise outer boundary conditions . 
This is demonstrated for instance in Fig. 5 of lNaravan et al.ld 19971) . 
where three very different outer boundary conditions on the gas ro- 
tation and temperature all give pretty much identical solutions in 
the mi ddle zone. The same is true also for Be (Fig. 7 of the same 
paper). lYuan et al ] d2012al) have carried out hydrodynamic simu- 
lations of ADAFs where they find that Be of the accreting gas is 
mainly set by the outer boundary condition. It is possible that their 
models do not extend over a large enough range of radius to sample 
the self- similar zone. 

All the results presented here refer to a non-spinning BH. 
This is the simplest version of the ADAF problem, where there 
is no additional complication from central energy injection by a 
spinning BH. It is also the case that relates most directly to the- 
oretical work as well as to non-relativistic MHD simulations. In 
the case of ADAFs around spinning BHs, although a large frac- 
tion of the energy from the BH seems to go out in a relativistic jet 
dTchekhovskov et aDl201ll) . some of it presumably propagates into 
the accreting flow. This energy very likely wi ll induce extra mass 
loss, a s see n in the simulations des cribed by iTchekhovskov et aD 
(1201 lh and iMcKinnev et all 120121) . Sorting out the BH spin ef- 
fect from the intrinsic effect due to ADAF physics is left for future 
work. 

In addition to outflows, we have also described in this paper a 
preliminary analysis of convection. In brief, the ADAF/SANE sim- 
ulation shows no evidence of convective instability (Fig. [XT}, while 

10 In the low-M RIAF branch of ADAFs, it is believed that the gas 
is two-temperature with non -relativistic ions and relativistic electrons 
dNaravan & McClintockll2008l) . If we take T e /T* = 0.1, areasonable value 
for an ion-dominated ADAF, then we expect F = 1.61. In the simulations 
presented here we have set T = 5/3, which is c lose enough, although 
techni cally in the "unphysical region" discussed by Mignone & McKinney 
(2007). In the ADAF literature, T = 1.5 is often used, but this is because 
those models wish to include the effect of a tangled magnetic field, which 
has an effective T = 4/3. In numerical MHD simulations, the magnetic 
field is treated as an independent component, so we are only concerned 
with the gas. Any choice T > 1.6 is probably reasonable. 



the ADAF/MAD simulation is apparently unstable by the Hoiland 
criteria over a part of its steady state region (Fig. QjO. However, 
there is little evidence in the MAD simulation for actual turbulent 
convection. Hence we speculate that the ADAF/ MAD simulatio n 
is probably in a state of frustrated convection (|Pen et al.l 120031) . 
Based on our current results, we are inclined to think that convec- 
tion is unimportant in ADAFs, whether SANE or MAD, but this 
issue needs to be investigated in greater detail before one can be 
certain. In particular, it is important to sort out the effect of the 
magnetic stress, which is ignored in the Hoiland criteria. Also, it is 
possible that the acc retion flow is descr i bed b y something like the 
global ID models in I Abramowicz et aD (120021) . where the flow be- 
haves like a basic ADAF (no outflow, no convection) until a radius 
r ~ 35 ru = 70, and then switches to a CDAF. We do not have 
enough dynamic range in our ADAF/SANE simulation to rule this 
out. 

We note that there are some observationa l indications against 
strong mass loss in ADAFs. [Allen et al showed that a num- 

ber of low-luminosity AGN have radio jets with implied powers 
that are a reasonable fraction of accretion energy at the Bondi rate 
from t he surrounding interstellar medium. In fact, iMcNamara et al.l 
(201 l|) identified systems with Pj e t > MeondiC 2 , and argued that 
these jets must be powered by BH spin. While it is true that a 
rapidly spinning BH can produce a very strong jet, the jet power 
is still linked to the accretion power; Pj et may be a factor of a 
few l a rger than Mbhc 2 , but not much more dTchekhovskov et all 
12011 : McKinn ev et aD 120121) . Therefore, the observations men- 
tioned above mean that a good fract ion of the available mas s at 
the Bondi radius must reach the BH dNaravan & Fabianll201ll) . If 
mass loss between the Bondi radius and the BH is very large, as in 
some versions of the ADIOS model ( Blandford & Begelman 19991 : 
lBegelmaDl2012h . or if a CDAF is present over a wide range of 
radius, there would not be sufficient mass near the BH to tap the 
BH spin energy and power the observed jets. We believe that the 
above observational evidence, assuming it holds up, drives us to- 
wards one of the following descriptions of the accretion flow: (i) 
an ADAF with a weak outflow, i.e., a value of the index s close 
to 0, or (ii) an ADAF with a strong outflow (s « 1) but with the 
outflow restricted to a small range of radius, say no more than one 
or two decades, or (iii) a CDAF with properties and scalings rather 
different from the analytical models in the literature dNaravan et al.l 
2000; Ouataert & Gruzinov 2000b), or (iv) a perfectly spherically 
symmetric Bondi flow. We consider the fourth possibility unlikely 
since it requires gas at the Bondi radius to have an extremely low 
specific angular momentum. 

The interesting differences we find between the ADAF/SANE 
and ADAF/MAD simulations brings up the question of which is 
more relevant for real systems. The defining feature of a MAD 
system is that accretion has dragged in a considerable amount of 
magnetic flux and has caused the field to accumulate around the 
BH. Whether or not ac cretion can drag field so effectively has 
been much debated (e.g.jLovelace. Rothstein & Bisnovatvi-Koganl 
120091 : iGuilet & Ogilvid l20ll and references therein), but it is 
agreed that field-dragging will be most efficient in thick accre- 
tion flows such as ADAFs rather than in thin discs. Assuming 
that inward advection of magnetic field does operate effectively in 
ADAFs, there is typically more than enough magnetic field avail- 
able i n the external mediu m to drive an accreting BH to the MAD 
state dNaravan et al.ll2003l) 
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